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The stationary functional of the all-electron density functional plus dynamical mean field theory 
(DFT+DMFT) formalism to perform free energy calculations and structural relaxations is imple¬ 
mented for the first time. Here, the first order error in the density leads to a much smaller, second 
order error in the free energy. The method is applied to several well known correlated materials; 
metallic SrV 03 , Mott insulating FeO, and elemental Cerium, to show that it predicts the lattice 
constants with very high accuracy. In Cerium, we show that our method predicts the iso-structural 
transition between the a and 7 phases, and resolve the long standing controversy in the driving 
mechanism of this transit on. 

PACS numbers: 71.27.+a,71.30.+h 


Prediction of the crystal structures of solids by large 
scale quantum mechanical simulations is one of the fun¬ 
damental problems of condensed matter physics, and oc¬ 
cupies a central place in materials design. The workhorse 
of the field is the Density Functional Theory (DFT) pQ at 
the level of Local Density Approximation (LDA) or Gen¬ 
eralized Gradient Approximations (GGAs), which pre¬ 
dict lattice constants of weakly correlated materials typ¬ 
ically within ^1% relative error [2|. 

These errors of DFT in LDA/GGA implementations 
are an order of magnitude larger in the so called cor¬ 
related materials: For example, the lattice constant of 
S -Pu is underestimated by 11% |3j or non-magnetic FeO 
by 7%[I]. While GGAs and hybrid functionals can some¬ 
times improve upon conventional LDA, these function¬ 
als many times degrade the agreement between predicted 
and experimentally determined bulk moduli and lattice 
constants, in particular in materials containing heavy el¬ 
ements. [2] 

To account for the correlation effects, more sophisti¬ 
cated many body methods have been developed. Among 
them, one of the most successful algorithms is the dy¬ 
namical mean-field theory (DMFT) [5]. It replaces the 
problem of describing correlation effects in a periodic 
lattice by a strongly interacting impurity coupled to a 
self-consistent bath |6y To become material specific, 
DMFT was soon developed into an electronic structure 
tool (LDA+DMFT) 00, which achieved great success 
in numerous correlated materials (for a review see [9]). 
The LDA+DMFT method has mainly been used for the 
calculation of spectroscopic quantities, and only a few 
dozens fT0M30] of studies managed to compute energetics 
of correlated solids, and only a handful of them used ex¬ 
act solvers and charge self-consistency na na eh EH Eg 
[293. This is not only because of the very high computa¬ 
tional cost, but also because previous implementations of 
LDA+DMFT were not stationary, and hence it was hard 
to achieve precision of free energies needed for structure 
optimization and study of phase transitions in solids. 

Here we implemented the LDA+DMFT functional, 


which delivers stationary free energies at finite temper¬ 
atures. This stationarity is crucial for practical imple¬ 
mentation and precision of computed energies, since the 
first order error in the density p (or the Green’s function) 
leads only to the much smaller second order error in the 
free energy, since the first order variation vanishes, i.e., 
SF/Sp = 0. This property is also crucial in calculating 
the forces, as stationarity of the functional ensures that 
only Hellmann-Feynman forces need to be computed for 
structural relaxation . [58] 

The DFT+DMFT total energy is given by [9] : 

E = Tr{H 0 G) + I'fr(EG) + E H [p] + E xc [p} 

^ |+Zoc] T Enuc—nuc (1) 

where i7 0 = —V 2 + £(r — r')V ext (r ), G is the elec¬ 
tron Green’s function, E H [p\ and E xc [p\ are Hartree 
and DFT exchange-correlation functional, V ext is the 
electron-nuclear potential, E nuc _ nuc is the interaction en¬ 
ergy of nuclei, E is the DMFT self-energy, and & DC [ni oc \ 
is the double-counting (DC) functional. [4] Here the 
Migdal-Galitskii formula (MGF) is used E pot = ^Tr(EG) 
to compute the DMFT part of the potential energy. 

Gordon Baym showed [31] that for certain class of ap¬ 
proximations, which are derivable from a functional ex¬ 
pressed in terms of closed-loop Feynman diagrams, MGF 
can be used instead of more complicated expression for 
evaluating the Luttinger-Ward Functional [32] [33] . He 
called such approximations conserving. While the DMFT 
is a conserving approximation in Baym’s sense, LDA or 
GGA are not, as the Galitskii-Migdal formula ^Tr (V xc p) 
has to be replaced by the exchange-correlation functional 
E xc [p\. As a result, the combination of DFT+DMFT in 
its charge-self consistent version is not conserving either, 
and consequently MGF can give different total energy 
than the Luttinger-Ward functional. Only the evaluation 
of the latter is guaranteed to give stationary free ener¬ 
gies. We will give numerical evidence that evaluation of 
MGF in Eq. [l] gives different results than evaluation of 
the Luttinger-Ward functional, which strongly suggests 
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that Eq.[l] gives non-stationary total energies. 

The Luttinger-Ward functional of DFT+DMFT has 
been well known for several years [9], but it has never 
been successfully implemented to compute the free en¬ 
ergy of a solids. It has the following form 


where Z atorn can be directly computed from atomic en¬ 
ergies, and Po is the probability for no kinks on the im¬ 
purity. Finally, Z = exp (—Fi mp /T), giving directly 

Fimp -T(log(Z atom ) - log(Po)). (5) 


r[G] = TrlogG - Tr((Gg 1 - G - 1 )G) + E H [p] 
+E xc [p\ + <Z> dmft [PG} - <S> DC [Pp\ + E nuc . nuc , ( 2 ) 

where G+ (it'; iui) = [iuj + p + V 2 - V ex t(r)]5(r - 
r'), $ DMFT [PG] is the DMFT functional, which is 
the sum of all local skeleton Feynman diagrams. 
The projected Green’s function PG = Gi oca i = 
T,lu\^l)^l\G\H'){H'\ and the projected density 
Pp = pi oca i are computed with projection to a set of lo¬ 
calized functions | 0 ) centered on the ’’correlated” atom. 
The projection defines the local Green’s function Gi OC ah 
the essential variable of the DMFT. 

The variation of functional T[G] with respect to G 
(5T[G]/5G) gives, 

G - 1 - Gq 1 + (V H + V xc )6 (r - r')5(r - r') 

y^ Mfr [gwi 

SG local 

_p 8$ DC \Plocal] S ( r _ r , )(5(r _ T ,) = Q) ( 3 ) 

0 Piocai 

which vanishes, since it is equal to the Dyson equation 
that determines self-consistent G, hence the functional is 
stationary. 

The value of the functional T at the self-consistently 
determined G delivers the free energy of the system M- 
We evaluate it by inserting G 0 1 — G 1 from Eq. [ 3 ] into 
Eq. [2] to obtain 

F = E nuc _ nuc - TY((Vh + V xc )p) + E H \p\ + E xc \p\ 

+Tr log G — Tr log G[ oc + Fi mp 
+Tr(V dc p loc ) - <S> DC [p loc ] + pN, (4) 

where we denoted V dc = 5$ DC [pi oca i\ /<5 pi oca i and Fimp 
is the free energy of the impurity problem, i.e., Fimp = 
Trlog Gioc - Tr(£G Zoc ) + $> DMFT [G loc \. [4] Here we also 
use the fact the solution of the auxiliary impurity prob¬ 
lem delivers the exact local Green’s function, i.e., E = 
5& DMFT [Glocal]/SGlocal, and we added pN because we 
work at constant electron number. 

The crucial point is that the continuous time quantum 
Monte Carlo method (CTQMC) [341135] solves the quan¬ 
tum impurity model (QIM) numerically exactly, hence, 
we can compute very precisely the impurity internal en¬ 
ergy as well as the free energy Fi mp of this model. At high 
enough temperature, Fi mp can be directly read-off from 
the probability for the perturbation order k , which we 
call Pk , and is computed by P^ = Z^jZ (where Z^ is the 
partition function with k— kinks), hence Z = Z a tom/F, o, 


When the temperature is low, Po becomes exponen¬ 
tially small, and we can no longer determine Z to high 
enough precision in this way. However, we can compute 
very precisely the internal energy of the impurity at ar¬ 
bitrary temperature. The internal energy of QIM E irnp 
is given by 


which follows directly from the thermodynamic av¬ 
erage of QIM Hamiltonian. Here the hybridization 
A and impurity levels Simp are determined from the 
local green’s function by the standard DMFT self- 
consistency condition G z “J aZ = iuj n — £i mp — £ — A, 
and Eimp-pot = |Tr (EGimp)- These quantities can be 
computed very precisely by CTQMC using the follow¬ 
ing tricks: i) Tr(AG^ mp ) is computed from the average 
perturbation order (k) of CTQMC, and takes the form 
Tr(AGi mp ) = (k) /T, where T is temperature [34]; ii) 
Fimp-pot is computed from the energies of atomic state 
of QIM E^ om and their probabilities P m by E irn p_p 0 t = 
Y,m P mEm° m - Tr (SimpTiimp) [Mj, which delivers much 
more precise interaction energy than obtained by MGF; 
ii) We spline A (cj n ) in Matsubara points and determine 
its derivative dA/doj n , and then carry out Matsubara 
sum by subtracting out the leading high-frequency tails 
by formula A/((iuj — £i)(iuj — £ 2 )), which has an analytic 
sum of A(f(e 1 ) — f{s 2 ))/(m — £ 2 )- Because probabilities 
P m and perturbation order (k) are known to very high 
precision in CTQMC, the impurity internal energy can 
easily be computed with precision of a fraction of a meV. 

To compute precise impurity free energy Fi mp at lower 
temperature, we first converge DFT+DMFT equations 
to high accuracy at low temperature. Using converged 
impurity hybridization A (icj n ), we raise the temperature 
of the impurity (keeping A fixed) to T>, so that Po is 
of the order of 10 -5 or higher, and obtain reliable Fimp 
using Eq. [5j and entropy Simp at this higher tempera¬ 
ture S> = (Eimp (T>) — Fimp(Ty ))/T>. Next, we evalute 
impurity internal energy for several inverse temperatures 
/3 = 1/T, and than we use standard thermodynamic re¬ 
lations to obtain entropy at lower temperature T by 


S(T) = S> - 


Fimp(F> ) Eimp(T) f 1/T 


nL/l 

J 1/T> 


Eimp(f3)d/3(7) 


where /3 = 1/T. This formula is obtained integrating 
by parts the standard formula S = f c v /TdT and c v = 
dE/dT. We hence obtain Simp and Fi mp = Ei m p — TSimp 
at lower T which can be inserted into Eq. [4] The rest 
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of the terms in Eq. [4] are relatively straightforward to 
evaluate, however, for a high precision implementation 
one needs to combine the terms that largely cancel and 
evaluate them together |4] . 

Previous implementations of free energy within 
LDA+DMFT [[25] [27] [36] were based on i) evaluating 
the total energy Eq. [l]at range of temperatures, and in¬ 
tegrating resulting specific heat [36], and ii) the coupling 
constant integration [251127 ], where total energy of the 
solid is needed for a range of coulomb repulsion’s U and 
is than integrated over U. In both approaches, the self- 
consistent LDA+DMFT solution is needed for many val¬ 
ues of the parameters (either U or T) to evaluate F. In 
our method, a single LDA+DMFT calculation for solid 
is needed, which makes the method much more efficient. 
Furthermore, current implementation of the free energy 
is stationary, hence higher precision of F is achieved. 

To test the implementation of the LDA+DMFT func¬ 
tional, we computed the volume dependence of the free 
energy for three well studied correlated materials: a 
metallic early transition metal oxide with perovskite 
structure SrVOs, a Mott insulating transition metal ox¬ 
ide FeO in its rock salt structure, and the lanthanide 
elemental metal, Cerium, in its face centered cubic struc¬ 
ture, which undergoes a first order iso-structural transi¬ 
tion. 

We used the implementation of LDA+DMFT of 
Ref|37, which is based on the Wien2K package [38] , and 
LDA in combination with nominal double-counting [371 
139] [40]. More technical details are given in the supple¬ 
mentary information. 

In Fig. [lja) we show the energy E(V), and F(V) 
for SrV 03 at T = 230 K, computed with Eq. |TJ and 
Eq. [4} respectively. The minima of E(V) and F(E) 
are achieved at 55.71 A and 55.51 A . The experimen¬ 
tally determined volume is V exp = 56.53 A [42]. The 
LDA+DMFT hence slightly underestimates the equilib¬ 
rium volume (1.8%), which gives 0.6% error in lattice 
constant. This is well within the standard error of best 
DFT functionals for weakly correlated materials. 

The metallic nature of SrVOs, with moderate mass en¬ 
hancements m* / mb anc i ~ 2 [4], leads to very small DMFT 
corrections in crystal structure [4|. Note that energy min¬ 
imization leads to slightly larger volume than free energy 
minimization, contrary to expectations. This is because 
energy is computed from non-stationary Eq. [l] while free- 
energy is obtained from the stationary expression Eq. [4] 
The latter is hence more trustworthy, and should be 
considered best LDA+DMFT result. This is also clear 
from pressure versus volume diagram in Fig. [IJb), where 
— dF/dV agrees more favorably with the experiment than 
- dE/dV obtained by MGF. 

In Fig. Qc), we show the impurity entropy obtained 
by Eq. 0 for two representative volumes. In this itin¬ 
erant system with very large hybridization, we do not 


Volume[A 3 ] 
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FIG. 1: (Color online): a) E(V) and F(V) for SrVOs at 
T — 230 K from Eq. [l] and [4] respectively. Entropy term 
TSimpiV ) is very small, (b) theoretical and experimental [41] 
p(V). Good agreement between theoretical —dF/dV and ex¬ 
periment is found, (c) Impurity Entropy Eq. [7] for represen¬ 
tative volumes. To obtain Simp, temperature is varied in the 
impurity problem only, and not in the LDA+DMFT problem 
of the solid. 

notice a temperature scale at which t2g shell is degen¬ 
erate (log (6)) nor the scale of the lowest order Kramers 
doublet (log( 2)), but we notice the Fermi liquid scale in 
the steep downturn of S(T) at T « 350 K. 

Fig.[2][a) shows E(V) and F(V) for paramagnetic Mott 
insulating FeO at 300 K, above its antiferromagnetic or¬ 
dering temperature. The equilibrium volume of E and 
F is 20.28A 3 and 20.24A 3 , while the experimental vol¬ 
ume is 20.342A . The lattice constant is thus under¬ 
estimated for only 0.10% and 0.16%, respectively. In 
comparison, all standard DFT functionals severally un¬ 
derestimate FeO lattice constant, for example PBE-sol, 
PBE, and LDA for 5.2%, 5.0%, and 7.7%, respectively. 

In Fig. [2|b) we show P(V) diagram and its excellent 
agreement with experiment. Fig. |2][c) shows impurity 
entropy Si mp (T ) for a few volumes. In contrast to metal¬ 
lic SrVOs, here we clearly see an extended plateau of 
Simp(T) = log(6)*kB around 1000 K, which signals com¬ 
plete degeneracy of the t^g shell, and its slight decrease 
at 300 K in proximity to the AFM state. 

The iso-structural transitions of Cerium attracted a 
lot of experimental and theoretical effort, but its theo¬ 
retical understanding is still controversial. On the ba¬ 
sis of LDA+DMFT calculation McMahan et.al nn pro¬ 
posed that the total energy exhibits a double-minimum 
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FIG. 2: (Color online): a) E(V) and F(V) for FeO from 
Eq.0 and [4] respectively. Entropy term TSimpiV ) is large 
but almost constant, (b) theoretical and experimental p(V). 
Filled and empty circles are from Refs.l43land l441 respectively, 
(c) Impurity entropy Eq. [t| for representative volumes. The 
degeneracy of the t 2 g shell above 1000K is apparent. 



FIG. 3: (Color online): a) E(V) and F(V) for elemental 
Cerium from Eq.[l]and [4] respectively. Data are presented for 
T=400 and 900K. (b) Entropy Simp(V) is large and changes 
dramatically accros the transiton. (c) theoretical and experi¬ 
mental |l5] p(V) diagram. 


shape, concomitant with the appearance of the quasipar¬ 
ticle peak at temperature as high as 1500 K, signaling the 
first order transition. Using different implementation of 
the same method, Amadon et.al 1271116] proposed that 
the transition is entropy driven, and that the total en¬ 
ergy is featureless with the minimum corresponding to 
low volume <a-phase. Only the addition of the entropy 
term moves the minimum to the larger volume of 7 -phase. 
In this picture the transition at low temperatures, where 
the entropy becomes small and cannot drive the tran¬ 
sition, is intrinsically absent. Yet another proposal was 
recently put forward on the basis of LDA+Gutzwiller cal¬ 
culations mm, in which the transition is present even 
at zero temperature, but the transition occurs at negative 
pressure. The transition is thus detectable even in the to¬ 
tal energy, in the absence of entropy, and becomes second 
order at T = 0. In the same method, the finite temper¬ 
ature transition is first order, and the double-minimum 
shape of free energy becomes most pronounced at very 
high temperature (1500 K) [48] . 

Our LDA+DMFT results for Ce are plotted in Fig. [3] 
The total energy curve at 400 K clearly shows a region 
of very flat shape in the region between the a-j volume. 
Indeed the derivative of the energy —dE/dV displayed 
in Fig. [3jc) shows a clear region of zero slope around 
1 GPa. This is consistent with results of Lanata et al. m 
finding very similar zero slope of —dE/dV at zero tem¬ 
perature, but is inconsistent with Ref. 27. which finds 
no feature in total energy. It is also inconsistent with 
McMahan et.al m showing clear double-peak in total 
energy. On the other hand, the addition of entropy sub¬ 
stantially increase the region of soft volume, as suggested 
by Amadon et.al [46] . Indeed the change of the entropy 
between the two phases is of the order of 0.9which 
is consistent with experimental estimations of 30 meV at 
400K [49]. The physical mechanism behind this large 
entropy change and unusual volume dependence of en¬ 
ergy is in very fast variation of coherence temperature, 
as suggested in Refs. P31I16], and conjectured in Kondo 
volume collapse theory m- The phase transition in our 
calculation occurs around 1.6 GPa, which is not far from 
experimentally determined critical pressure of 1.25 GPa 
at T = 400 K. The free energy barrier in our calculation 
is however extremely small, and no clear double peak 
of F(V) or negative slope of —dF/dV can be detected 
within our 1 meV precision of energies. This is similar to 
results of Ref. 08] at 400 K, but different from Ref. [TT] 
While the start of the transition region in <a-phase is 
in good agreement with experiment, the 7 -phase vol¬ 
ume is underestimated in our calculation. We believe 
that the addition of phonon entropy is needed to further 
increase the transition region, and establish larger free 
energy barrier between the two phases. Experimentally, 
above 460 K the a — 7 phase transition ends with the fi¬ 
nite temperature critical point. Our calculation at high 
temperature 900 K shows that the signature of the phase 
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transition in F(V) and E(V) disappears, which is dif¬ 
ferent than predicted by Gutzwiller method [48], where 
the largest free energy barrier is found at these elevated 
temperatures, but qualitatively consistent with Ref. QT| . 

In summary, we successfully implemented the station¬ 
ary formula for the free energy of DFT+DMFT method. 
On the example of SrVOs, FeO and Ce metal we demon¬ 
strated that the method successfully predicts lattice vol¬ 
umes in correlated solids, which are difficult for standard 
DFT functionals. We also resolved controversy in the 
mechanism of the a -7 transition in Cerium. 
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TECHNICAL DESCRIPTION OF 
COMPUTATIONAL DETAILS 

We used the implementation of LDA+DMFT of 
RefE2 which is based on Wien2K package [38] . The 
exchange-correlation functional E xc of LDA was utilized, 
in combination with nominal double-counting (DC) [37l 
139] . which was shown to be closest to the exact form of 
DC [3D]. We checked (on the example of Cerium) that the 
exact-DC gives very similar free energy, as expected for a 
stationary functional. The convergence of LDA+DMFT 
results is much faster using nominal DC, hence most of 
results in this publication are obtained by this simplifi¬ 
cation. 

The impurity model is solved using the hybridization 
expansion version of the numerically exact continuous 
time QMC method [34] 35]. Of the order of 300 LDA 
and 30 DMFT iterations were required for precision of 
1 meV per formula unit, and between 10 9 — 10 10 Monte 
Carlo moves were accepted per impurity iteration for pre¬ 
cise enough impurity solution. The resources of Titan 
supercomputer were used. 

To construct the projector, the atomic-like local or¬ 
bitals are used (r|= ^^V m (r). The radial part of 
the local orbital ui(r) is the solution of the scalar rela¬ 
tivistic Dirac equation inside the muffin-tin sphere, lin¬ 
earized at the Fermi level. The muffin-tin spheres are 
set to touch at the lowest volume. We tested a few other 
forms of the projectors defined in Ref. 37. The stationary 
F(V) is quite insensitive to the precise choice of projec¬ 
tor, however, E(V) changes much more. 

For SrVOs calculations, we treated dynamically all five 


3d orbitals of Vanadium. The muffin-thin radius of Vana¬ 
dium was set to R rnt = 1.83a#, and U = 10eV was 
used, which was previously shown to give good spectra 
[39] [40] for this localized orbital (see spectra below). The 
Yukawa form of screening interaction than gives J ~ 1 eV 
(see note below). Brillouin zone integrations were done 
over 15x15x15 k-point in the whole zone in the self- 
consistent calculations, and for calculation of the impu¬ 
rity entropy, the hybridization is computed on more pre¬ 
cise 36 x 36 x 36 k-points mesh. We mention in pass¬ 
ing that impurity entropy is very sensitive to the precise 
frequency dependence of the hybridization, and requires 
very dense momentum mesh. 

For FeO, all five 3d orbitals are treated by DMFT and 
the muffin-thin radius of iron is set to R mt = 2.11 a#, 
and the Coulomb repulsion to previously determined V = 
SeV [5F, which requires J ~ leV in Yukawa form. In 
Ce metal, all seven 4/ orbitals are treated by DMFT and 
the muffin-thin sphere is Rmt = 2.5 a#, the k-point mesh 
is 21 x 21 x 21, and the Coulomb V = 6eV [27l 1461152] . 
leads to J = 0.72 eV in Yukawa form. The spin-orbit 
coupling is included in Cerium, but neglected in SrVOs 
and FeO. 


DETAILS ON EVALUATION OF LDA+DMFT 
FUNCTIONAL 

Here we explain how we evaluate the total energy Eq. 1 
and the free energy Eq.4 of the main text. 

For total energy Eq.l, we group the terms in the fol¬ 
lowing way 


E = TY((-V 2 + V ext + V H + V XC )G) - TV((Vy + V xc )p) + E H [p\ + E xc [p] + E nuc . nuc + ^Tr(EG) - ^> DC [pi oc ] (8) 


We then split the energy into three terms E = E\ + E^ + 
F? 3 , where the first two parts Ed, E 2 are computed using 
the Green’s function of the solid, and the third E% using 
the impurity Green’s function. 

The first five terms in Eq.[8]look similar to the standard 
DFT energy functional, except that the Green’s func¬ 
tion G here is the self-consistent LDA+DMFT Green’s 
function. We first solve the eigenvalue problem for 
Kohn-Sham states (—V 2 + V ex t + Vh + V xc )^u c = 
where ef k FT are DFT-like energies, computed 


on LDA+DMFT charge. We then evaluate 

Ei = Tt(£ DFT G) (9) 

and 

E 2 = -Ti((Vh + v xc )p ) + E H [p\ + E xc [p] + E nu _ nu ll0) 

Both Ei and E 2 are computed using Green’s function G 
and density p of the solid in the same way as the standard 
DFT total energy is implemented [53] . 

The last two terms of Eq. [8] can be computed either 
from the local Green’s function PG or from the impu- 





rity Green’s function Gimp. Once the self-consistency is 
reached, the two are of course equal. We choose to eval¬ 
uate the second term on the impurity Gimp 

F 3 = ^^(J^impGimp) — & DC [Pimp] (H) 

However, we never actually use Migdal-Galitskii formula, The free energy functional T[G] (Eq. 2 of the main 
because it is numerically much less stable than computing text) is 


the potential energy from the impurity probabilities, i.e., 
^Tr(I \mpGimp) = ^PmEm 0 ™ - Tr (e imp n imp ) 


r[G] = Tr log G — lVlog((Go 1 - G~ 1 )G) + E H [p ] + E xc \p] + <!> DMFT [G loc ] - $ DC [p loc ] + E nuc _ nuc . (12) 


First, we extremize it (ST[G]/SG = 0) to obtain the 
Dyson equation 


G — G 0 + Vh + V xc + S dm ft ~ V dc = 0. (13) 

A note is in order here. We assumed SP/SG = 0, which 
holds whenever the projector does not depend on the self- 
consistent charge density. To ensure this property, we 
used for the localized orbitals \<p) = -^-^T/ m (f), where 
the radial wave function ui(r ) is the solution of the scalar 
relativistic Dirac equation on the LDA charge density 
(rather than self-consistent charge density). Note also 
that the use of the self-consistently determined Wannier 
functions (which depend on self-consistent charge), as is 
commonly used in most of the LDA+DMFT implemen¬ 
tations [ 2U E3 EQ] , leads to non-stationary LDA+DMFT 
solution, and non-stationary free energies. 

We next insert the expression G _1 
to obtain expression for free energy 


G 0 1 into Eq. 


12 


F = E nuc — nuc - Tt((V h + V xc )p) + E H \p\ + E xc [p\ 

+Tr log G — Tr(X DMFT G) + <S> DMFT [G loc \ 

+Tr(V dc pioc) ~ $ DC [pioc} + pN(14) 

The impurity free energy Fi mp contains & DMFT [Gi m p\ in 
the following way 


Fimp Tr log Gimp Fl(EiimpG im p 


)+$> DMFT [G l mp[l 5 ) 


In DMFT, Gioc — Gimp nnd — Fjimpi hence we 

can write 


F = E nuc _ nuc - Tr((V H + V xc )p) + E H \p\ + E xc [p] 
+Tr log(G) - Tr log(G; oc ) + F imp 
+Tr (Vdcpioc) ~ & DC [pioc\ + //.V.( l(i) 

This equation appears as Eq.4 in the main text. 

Next we split free energy of the impurity into the en¬ 
ergy and the entropy term 


Fimp — Fimp FSimpi 


where 


Fimp IV((A + Simp dJn t )Gimp) 

dOJ n 

F-FT(TiimpG imp ) FSimp 


(17) 


Hence 


1 . 


F = -Tr (ZimpGimp) - [pioc] - TSi % 


+E nuc _ nuc - Tr((V H + V xc )p) + E H [p\ + E xc [p\ 

+Tr log(G) - Tr log(G loc ) + Tr(V dcPloc ) + pN 

+Tr((A T £imp ^ n ~j (1$) 

cL(jj n 

Again using the identity G^p = Gi oc and p im p = Pioc as 
well as the definition of E 2 (Eq. [To]) and E3 (Eq. 11) we 
obtain 


F — Tr log(G) T pN T E 2 

T Tr((A uj n — (- Simp T Vdc)Gioc) 

duj n 

— Tr log (Gioc) + E 3 — TSimp (19) 

This equation is implemented in our DFT-fDMFT code. 
Similarly than in the implementation of the total energy 
Eq. [8j we compute E 3 and TSimp using impurity quanti¬ 
ties, while the rest of the terms are computed using the 
Green’s function of the solid. [Alternatively, we could 
also compute the last two rows using impurity quantities, 
and the first row using the solid Green’s function]. In this 
way we ensure that F and E are split in the same way be¬ 
tween the ’’impurity” and the ’’lattice” quantities, hence 
they share almost identical Monte Carlo noise. However, 
when comparing E(V) and F(V) at two different vol¬ 
umes, F(V) converges faster that E(V) with the number 
of LDA and/or DMFT iterations. Moreover, F(V) is 
very robust with respect to small changes in projector or 
double-counting, while E is more sensitive. 

Notice that F + TSi mp can be evaluated at each 
LDA+DMFT iteration, just like the total energy above. 
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To add TSimp at low temperatures, we however need 
a few extra impurity runs. The method of computing 
TSimp is explained in the main text, and requires the 
impurity energy at a few temperatures. An alternative 
to this approach is to compute TSimp from so called 
” flat-histogram sampling method” [54], which is also 
done as postprocessing on self-consistent LDA+DMFT 
hybridization A. 

Perhaps, the most challenging term in Eq. [l9]to com¬ 


pute is Trlog(G), which requires eigenvalues (but not 
eigenvectors) of the LDA+DMFT eigenvalue problem. 
We first diagonalize 

(-V 2 + V ext + Vh + V xc + Yi(iu n ) — Vdc)^i,k,u n = 

— £i,k,(jj n '*Pi,k,uj n • ( 20 ) 

and then evaluate 


Tr log(G) + iiN = T ^ (log(£ i)fc)W „ - - ft)- log(£; jfei00 - iu n 

iu n ,i,k,cr 


h))-tJ2 log(l + + nN( 21) 

i,k,a 



FIG. 4: Free energy of LDA+DMFT for SrVOs compared 
with total energy of other standard DFT functionals. 


Here it becomes apparent that if E(icj n ) is frequency in¬ 
dependent, the first term in the brackets vanishes, while 
the second term gives (at T = 0) the sum of eigenvalues 

Tr log(G) + fiN -+ u ~°—>• 0(£i,k < aO 

i,k,cr 

the well known DFT contribution to the total energy. 



FIG. 5: Free energy of LDA+DMFT for FeO compared 

with total energy of other standard DFT functionals. Upper 
(lower) panel shows non-magnetic (antiferromagnetic) DFT 
calculation. LDA+DMFT results are obtained at 300K in 
paramagnetic state. 


COMPARISON WITH STANDARD 
FUNCTIONALS 

Here we compare total energy of LDA, PBE [[55] . 
and PBEsol [56] functionals with the free energy of 
LDA+DMFT. 

In most weakly correlated solids, LDA underestimates 
lattice constants on average for 1.6%, while PBE [55] 
overestimates them for approximately 1%. [2] PBEsol [56] 
was designed to predict most accurate volumes in solids, 
and it typically falls in-between LDA and PBE. 

In Fig. [4] we compare LDA+DMFT free energy in 
SrVOs with the total energy computed by other function¬ 
als. Both LDA+DMFT and PBEsol underestimate lat¬ 
tice constant for approximately 0.6%, while LDA under¬ 


estimates it for 1.5%, and PBE overestimates for 0.7%. 
Hence predictions of standard functionals in the case of 
SrVOs are quite in line with standard performance in 
weakly correlated solids. Perhaps, this is not very sur¬ 
prising given that SrV 03 is a metallic moderately corre¬ 
lated system. 

In FeO (Fig. [ 5 ]), all standard functionals severally un¬ 
derestimate volume in the paramagnetic state. For ex¬ 
ample the lattice constants with LDA, PBEsol and PBE 
are 7.7%, 6.5% and 5.1% too small, far outside the stan¬ 
dard performance of these functionals in weakly corre¬ 
lated solids. 

The predictions are improved when the AFM long 
range order is allowed. LDA and PBEsol still underes¬ 
timate lattice constant for 3.6%, and 2.3% respectively. 
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FIG. 6: Free energy of LDA+DMFT for Cerium compared 
with total energy of other DFT functionals. LDA+DMFT 
results are obtained at 400 K. 

On the other hand PBE is this time quite close to the 
experiment (underestimates for 0.7%). In comparison 
LDA+DMFT underestimates it for only 0.16%. It is 
quite clear that the excellent prediction of AFM-PBE 
here is merely a coincidence, as normally PBE overesti¬ 
mates the volume. 

Finally, we plot results for Cerium in Fig. [6] The result 
of LDA+DMFT is very different from those of any other 
functional, as it clearly contains the nontrivial soft mode 
for the a-j transition. No other functional shows any 
hint of such transition. 


The equilibrium volume in Cerium is strongly temper¬ 
ature dependent, and is approximately 34A at zero pres¬ 
sure and 400 K, while it changes to approximately 28A 
in the a phase at low temperature. The LDA+DMFT 
results are computed at 400 K, hence at p = 0 the volume 
is somewhat underestimated (1.5%), but under pressure 
(already at 1 GPa) the agreement with experiment is con¬ 
siderably improved. 

The DFT results should be compared to T = 0 ex¬ 
perimental volume of 28 A . All functionals underesti¬ 
mate the lattice constant, LDA for 6%, PBEsol for 5% 
and PBE for 1.8%. Clearly electronic correlations are 
very important even in the a phase at low temperature, 
as standard DFT functionals substantially underestimate 
the volume. 


SCREENED COULOMB REPULSION OF 
YUKAWA FORM 

It is noted above that we used the Yukawa represen¬ 
tation of the screened Coulomb interaction, in which 
there is unique relationship between the Hubbard U and 
Hund’s coupling J. If U is specified, J is uniquely de¬ 
termined. To show this we derive the matrix elements of 
screened Coulomb interaction in our DMFT orbital basis 


U 7 




/+ 


,3 / ( M r )\ 2 ( ui(r') X 2 


d 6 r 


A|r rl 


r - r' 


There exist a well known expansion of Yukawa interaction in terms of spheric harmonics Y& m , which reads 

— A|r—r'l 

1 —V = 47r 

Ir-r'l ‘ 




( 22 ) 


(23) 


Here r< = ram(r, r'), r> = max(r, r'), I and K are modified Bessel function of the first and second kind. Inserting 
this expression into Eq. [22] we get 


47T 


^mim 2 m 3 m 4 — ^ ^ 2 k 1 ^ lrni I^Zm 4 ) 0 ^lm 2 II ^lm 3 ) 


pOO pOO 

:(2k + l) / dr dr'uf(r)uf(r') 

Jo Jo 


^fe+l/2(Ar<(Ar>) 


V r < r > 

Hence, the screened Coulomb interaction has the Slater form with the Slater integrals being 

/x A+l/2(^<)^£+l/2(^ r >) 


F _ 


P OO POO 

(2k + 1) / dr dr'uf(r)uf(r')- 

J 0 Jo 


V r < r > 


(24) 


(25) 


This is a product of two one-dimensional integrals and is 
very easy to efficiently implement. 


It is clear from Eq. [25] that A uniquely determines 
all F k ’ s, and furthermore even one Slater integral (F°) 
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FIG. 7: Spectral function of SrV 03 within LDA+DMFT at 
equilibrium volume compared with ARPES spectra of Ref.l57l 


uniquely determines A. This is because F k are monotonic 
functions of A and take the value of bare F k at A = 0 and 
vanish at large A. Hence given F°, the screening length 


A is uniquely determined, and hence other higher order 
F k are uniquely determined as well. 


MASS RENORMALIZATION OF METALLIC 

SrVOs 

Even though the Coulomb interaction in SrVOs is U = 
10 eV, it gives a relatively moderate mass enhancement 
over DFT band structure in all-electron LDA+DMFT 
implementation. This is because the interaction is 
severely screened by hybridization of d states with oxy¬ 
gen p states, and because the F g orbitals are in mixed- 
valence state (n t 2 g ~ 1.5) [39| [40]. In Fig. [7] we show the 
LDA+DMFT spectral function as well as recent APRES 
measurements m ■ The mass renormalization in the t 2 g 
orbital is ml 2g /m hand « 2 and in e g is m% 2g /m band « 1.3 
The agreement between ARPES spectra (the experimen¬ 
tal signal is color coded on the right) and LDA+DMFT 
spetrcal fuction A(fc, uj) (plotted on the left) is very good, 
both in the quasiparticle band (between —0.5eV and 
0.5 eV) and Hubbard sattelite at —1.5eV. 






